function [ DrM ] = getDrM( r, rn, prpn, M, qn )
%GETDRM Summary of this function goes here
%   Detailed explanation goes here

down = (r(1)^2+r(2)^2)^2;
up1 = r(1)^2-r(2)^2;
up2 = r(2)^2-r(1)^2;

s = -1 / rn * prpn;
A = 2*M;
ds = [(2*qn(2)*r(1)*r(2)+qn(1)*up1), (2*qn(1)*r(1)*r(2)+qn(2)*up2)] / down;
dA2 = [4*r(1)*r(2)^2, 2*r(2)*up2, -4*r(1)^2*r(2), 2*r(1)*up1; 2*r(2)*up2, -4*r(1)*r(2)^2, 2*r(1)*up1, 4*r(1)^2*r(2)] / down;
dA = reshape(dA2, 2, 2, 2);

DrM = tensorProductGradient(s, ds, A, dA);

end

